# average
avg_lambda = lambda_H * (lambda_L_to_H/(lambda_H_to_L+lambda_L_to_H)) + lambda_L * (lambda_H_to_L/(lambda_H_to_L+lambda_L_to_H))
7/212
kappa_lambda_fun <- function(lambda){
  kappa_lambda = (lambda-lambda_L)*(lambda_H-lambda)/lambda
  return(kappa_lambda)
}
mu_lambda_fun <- function(lambda){
  mu_lambda = (lambda_L-lambda)*lambda_H_to_L + (lambda_H-lambda)*lambda_L_to_H  - (lambda-lambda_L)*(lambda_H-lambda)
  return(mu_lambda)
}

# library(markovchain)
N_sim_years = 1000
dt = 1/12
simulation_lambda <- function(N_years=1000, dt=1/12, lambda_0=lambda_L, lambda_L_to_H=lambda_L_to_H){
  lambda = lambda_0
  lambda_sim = c(lambda)
  N_total = round(N_years/dt)
  random_vector = runif(N_total)
  for( i in 1:N_total ){
    if(lambda==lambda_L & (random_vector[i]<lambda_L_to_H*dt)){
      lambda = lambda_H 
    }else if(lambda==lambda_H & (random_vector[i]<lambda_H_to_L*dt)){
      lambda = lambda_L
    }
    lambda_sim = c(lambda_sim, lambda)
  }
  return(lambda_sim)
}

simulation_hat_lambda <- function(N_years=1000, dt=1/12, seed=1, lambda_L_to_H=0.1){
  TB_results = data.frame()
  N_total = N_years/dt
  if(seed!=-1){
    set.seed(seed)
    lambda_sim = simulation_lambda(N_years=N_years, dt=dt, lambda_0=lambda_L, lambda_L_to_H=lambda_L_to_H)
    dN_vec = runif(N_total)<lambda_sim*dt
  }else{
    lambda_sim = rep(0, N_total)
    dN_vec = c(rep(0, round(N_total/3)), 1, rep(0, round(N_total/3)-1), 1, rep(0, round(N_total/3)-1)  )
  }
  lambda = (lambda_H+lambda_L)/2
  init_vec = rep(0, N_total)
  TB_results = data.frame(i=init_vec, lambda_sim=init_vec, lambda=init_vec, dN=init_vec)
  for(i in 1:N_total){
    if( i%%100==0 ){
      print( paste0("Iteration ", i) )
    }
    kappa_lambda = kappa_lambda_fun(lambda)
    mu_lambda = mu_lambda_fun(lambda)
    dN = dN_vec[i]
    lambda = lambda + mu_lambda*dt + kappa_lambda*dN
    
    TB_results[i,] = c(i, lambda_sim[i], lambda, dN)
  }
  return(TB_results)
}


irrational_lambda <- function( lambda_vec, dt, lag, theta=1 ){
  # The unit of lag is in years
  # dt is the basic time unit in years
  N = round(lag/dt)
  a = lambda_H_to_L
  b = lambda_L_to_H
  
  aH = (b+a*exp(-(a+b)*lag)) / (a+b)
  bH = 1-aH
  aL = (b-b*exp(-(a+b)*lag)) / (a+b)
  bL = 1 - aL
  
  if(length(lambda_vec)<=2*N){
    print("Simulation length is too short!")
    return(NA)
  }else{
    lambda_vec_lag = c( rep(lambda_vec[1],N), lambda_vec[1:(length(lambda_vec)-N)] )
    lambdaT_lambdaL = aH*(lambda_vec_lag-lambda_L) + aL*(lambda_H-lambda_vec_lag)
    lambdaH_lambdaT = bH*(lambda_vec_lag-lambda_L) + bL*(lambda_H-lambda_vec_lag)
    ratio= (lambdaT_lambdaL/lambdaH_lambdaT) / ((lambda_vec-lambda_L)/(lambda_H-lambda_vec)) 
    lambda_theta_lambdaL = (lambda_vec-lambda_L)  *  (lambda_H-lambda_L)/( ratio^theta*(lambda_H-lambda_vec) + (lambda_vec-lambda_L) )
    lambda_theta = lambda_theta_lambdaL +  lambda_L
    return(lambda_theta)
  }
}





